Scenario

Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.

Load packages

library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement 
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed 
library(bayesplot) # needed for MCMC diagnostics 
library(DHARMa) # model validation 
library(ggdist) # partial plots 
library(tidybayes) # partial plots 
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans 
library(ggpubr) # arranging figures 
library(mgcv) # GAM 
library(modelr) # plotting

Set working directory

knitr::opts_knit$set(root.dir=working.dir)

Import data

This data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty

lat_resp_dat <- read_delim("C:/Users/jc527762/OneDrive - James Cook University/PhD dissertation/Data/Chapter4_Latitude/import_files/lat_resp_dat_no_outliers.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(area = col_skip(), rate.a.spec = col_skip()), 
    trim_ws = TRUE)

Data manipulation

rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |> 
  mutate(dev.temp = as.factor(dev.temp), 
         replicate = str_sub(sampleID, -1,-1), 
         population = factor(population)) |>
                                           # number of observations = 5758
  filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
         sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
         sampleID != "60.LCKM.152.30.1"  # 5637 - 76 = 5542
  ) |> 
  filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes) 
  group_by(sampleID) |> 
  mutate(max_value_index = which.max(rate.output), 
         row_number = row_number(), 
         max.rate = max(rate.output), 
         low.rate = mean(rate.output[order(rate.output)[1:10]]), 
         net.rate = max.rate - low.rate) |> 
  filter(row_number <= max_value_index) |>
  ungroup() |> 
  mutate(region = factor(region), 
         dev.temp =factor(dev.temp), 
         level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
                           tank >= 200 & tank <= 299 ~ 2,
                           tank >= 300 & tank <= 399 ~ 3,
                           TRUE ~ NA_real_)), 
         female = factor(female))

Exploratory data analysis

2nd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

3rd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

Fit model [random factors]

First we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.

Hypothesis test will include:

Null model
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK)
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL)
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)
EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)+ (1| CLUTCH_ORDER)

lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output)))

## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors <- prior(normal(430, 0.25), class="Intercept") +  
  prior(student_t(3, 0, 195), class = "sigma")

Models

Null

f.model.null <- bf(rate.output ~ 1, 
                   family = gaussian()) 

model.null <- brm(f.model.null,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model.null, file="./03.CTmax_resp_data_analsyis_files/model.null.RDS")

Model1

f.model1 <- bf(rate.output ~ 1 + (1| female) + (1| tank), 
                          family=gaussian())

model1 <- brm(f.model1,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1, file="./03.CTmax_resp_data_analsyis_files/model1.RDS")

Model2

f.model2 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level), 
                          family=gaussian())

model2 <- brm(f.model2,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 13 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model2, file="./03.CTmax_resp_data_analsyis_files/model2.RDS")

Model3

f.model3 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population), 
               family=gaussian())

model3 <- brm(f.model3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model3, file="./03.CTmax_resp_data_analsyis_files/model3.RDS")

Model4

f.model4 <- bf(rate.output ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order), 
                        family=gaussian())

model4 <- brm(f.model4,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

saveRDS(model4, file="./03.CTmax_resp_data_analsyis_files/model4.RDS")

LOO

LOO(model.null,model1, model2,model3,model4)

## Output of model 'model.null':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -31469.1  78.1
## p_loo         2.5   0.4
## looic     62938.2 156.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.6  99.6
## p_loo        57.6   3.1
## looic     61325.1 199.1
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.7  99.5
## p_loo        57.6   3.0
## looic     61325.4 198.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.7  99.6
## p_loo        58.1   3.2
## looic     61325.5 199.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model4':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.5  99.4
## p_loo        57.5   3.0
## looic     61325.0 198.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##            elpd_diff se_diff
## model4        0.0       0.0 
## model1       -0.1       0.4 
## model2       -0.2       0.4 
## model3       -0.2       0.4 
## model.null -806.6      48.2

Fit model [fixed factors - polynomials]

Prior investigation

lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output)))

## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors <- prior(normal(430, 0.25), class="Intercept") + 
  prior(normal(0, 40), class="b")
  prior(student_t(3, 0, 195), class = "sigma")

Models

Model1.p1

f.model1.p1 <- bf(rate.output ~ 1 + 
                         dev.temp*region + Temperature +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p1 <- brm(f.model1.p1,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.p1, file="./03.CTmax_resp_data_analsyis_files/model1.p1.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98596, p-value = 0.704
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p1$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p1$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p1$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p1$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p1$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

model1.p2

f.model1.p2 <- bf(rate.output ~ 1 + 
                         dev.temp*region + poly(Temperature, 2) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p2 <- brm(f.model1.p2,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.p2, file="./03.CTmax_resp_data_analsyis_files/model1.p2.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98971, p-value = 0.688
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p2$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p2$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p2$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p2$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p2$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

model1.p3

f.model1.p3 <- bf(rate.output ~ 1 + 
                         dev.temp*region + poly(Temperature, 3) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p3 <- brm(f.model1.p3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.p3, file="./03.CTmax_resp_data_analsyis_files/model1.p3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98512, p-value = 0.616
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p3$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

LOO

LOO(model1.p1, model1.p2,model1.p3)

## Output of model 'model1.p1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -29946.0  94.1
## p_loo        60.3   3.1
## looic     59892.1 188.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30525.1 100.7
## p_loo        58.9   3.2
## looic     61050.1 201.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30516.1 100.1
## p_loo        58.7   3.1
## looic     61032.1 200.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.2]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##           elpd_diff se_diff
## model1.p1    0.0       0.0 
## model1.p3 -570.0      31.3 
## model1.p2 -579.0      31.5

Fit model [fixed factors - GAM]

GAM model k=4

f.model1.gamk4 <- bf(rate.output ~ 1 + 
                         t2(dev.temp,region,Temperature, k=4, bs="fs") + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

model1.gamk4 <- brm(f.model1.gamk4,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.gamk4, file="./03.CTmax_resp_data_analsyis_files/model1.gamk4.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gamk4, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk4, ndraws=250, summary=FALSE)
model1.gamk4_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gamk4_resids) ; testDispersion(model1.gamk4_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97732, p-value = 0.464
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gamk4$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gamk4$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gamk4$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gamk4$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gamk4$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

### {-}

GAM model k=3

f.model1.gamk3 <- bf(rate.output ~ 1 + 
                         t2(dev.temp,region,Temperature, k=3, bs="fs") + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

model1.gamk3 <- brm(f.model1.gamk3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95))

## Compiling Stan program...

## Start sampling

saveRDS(model1.gamk3, file="./03.CTmax_resp_data_analsyis_files/model1.gamk3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gamk3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk3, ndraws=250, summary=FALSE)
model1.gamk3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gamk3_resids) ; testDispersion(model1.gamk3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9811, p-value = 0.536
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gamk3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gamk3$fit)

## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gamk3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gamk3$fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gamk3$fit, separate_chains = TRUE)

## 'pars' not specified. Showing first 10 parameters by default.

### {-}

AIC(model1.gamk3, model1.gamk4)

## Warning in AIC.default(model1.gamk3, model1.gamk4): models are not all fitted
## to the same number of observations

Partial effects plots

conditional_effects

model1.gamk3 |> 
  conditional_effects(spaghetti = TRUE, ndraws=200)

Model investigation

summary

summary(model1.gamk4)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rate.output ~ 1 + t2(dev.temp, region, Temperature, k = 4, bs = "fs") + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank) 
##    Data: lat_resp_dat2 (Number of observations: 4776) 
##   Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
##          total post-warmup draws = 1800
## 
## Smooth Terms: 
##                                    Estimate Est.Error l-95% CI u-95% CI Rhat
## sds(t2dev.tempregionTemperature_1)   809.91    168.24   551.16  1213.34 1.00
## sds(t2dev.tempregionTemperature_2)   878.19    186.51   593.72  1318.30 1.00
##                                    Bulk_ESS Tail_ESS
## sds(t2dev.tempregionTemperature_1)     1300     1719
## sds(t2dev.tempregionTemperature_2)     1681     1669
## 
## Group-Level Effects: 
## ~female (Number of levels: 15) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)   113.90     32.16    62.97   185.66 1.00     1535     1705
## 
## ~tank (Number of levels: 52) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)    82.15     10.73    64.41   107.17 1.00     1721     1643
## 
## Population-Level Effects: 
##                                  Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                          430.00      0.26   429.48   430.50 1.00
## scalemasscenterEQTRUEscaleEQTRUE    19.84      3.65    12.56    27.16 1.00
##                                  Bulk_ESS Tail_ESS
## Intercept                            1766     1635
## scalemasscenterEQTRUEscaleEQTRUE     1894     1643
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   114.30      1.18   111.91   116.61 1.00     1713     1567
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

#car::Anova(model1.geo.gamma)

R2

model1.gamk4 |> bayes_R2(summary = FALSE) |> median_hdci()

tidyMCMC

tidyMCMC(model1.gamk4, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

gather draws

#model1.re.wo |> get_variables()
model1.gamk4 |> gather_draws(`b_.*|sigma`, regex =TRUE) |> 
  median_hdci()

bayesplot

model1.gamk4 |> mcmc_plot(type='intervals')

emmeans - pariwise

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> regrid()

##  region  dev.temp emmean lower.HPD upper.HPD
##  core    28          407       325       502
##  leading 28          404       297       516
##  core    30          447       367       520
##  leading 30          387       268       494
##  core    31          452       360       531
##  leading 31          468       374       580
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()

emtrenns - pariwise

emtrends(model1.gamk4, specs=pairwise~dev.temp*region, var = "Temperature")

## $emtrends
##  dev.temp region  Temperature.trend lower.HPD upper.HPD
##  28       core                42.04     29.76      56.2
##  30       core                 1.39    -10.50      12.5
##  31       core                20.20      8.77      31.8
##  28       leading             -1.66    -15.93      11.0
##  30       leading              1.03    -13.20      14.5
##  31       leading              4.25     -8.60      16.4
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast                                estimate lower.HPD upper.HPD
##  dev.temp28 core - dev.temp30 core         40.607    24.078     60.12
##  dev.temp28 core - dev.temp31 core         21.926     3.890     39.59
##  dev.temp28 core - dev.temp28 leading      44.028    24.126     61.31
##  dev.temp28 core - dev.temp30 leading      41.138    21.869     60.16
##  dev.temp28 core - dev.temp31 leading      37.628    20.663     57.08
##  dev.temp30 core - dev.temp31 core        -18.947   -35.413     -2.84
##  dev.temp30 core - dev.temp28 leading       2.996   -14.071     20.71
##  dev.temp30 core - dev.temp30 leading       0.223   -17.503     19.71
##  dev.temp30 core - dev.temp31 leading      -3.028   -19.615     12.87
##  dev.temp31 core - dev.temp28 leading      21.958     4.079     38.88
##  dev.temp31 core - dev.temp30 leading      19.222    -0.815     36.84
##  dev.temp31 core - dev.temp31 leading      15.969    -1.138     32.29
##  dev.temp28 leading - dev.temp30 leading   -2.709   -23.119     15.26
##  dev.temp28 leading - dev.temp31 leading   -5.889   -23.118     13.74
##  dev.temp30 leading - dev.temp31 leading   -3.028   -20.894     16.54
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

probabilities

mtsqst <- model1.gamk4 |> emmeans(pairwise ~ region*dev.temp)
mtsqrt2 <- mtsqst$contrasts |> gather_emmeans_draws()
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())

summary figure

resp.plot <- lat_resp_dat2 |> 
  group_by(region, dev.temp) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~dev.temp); resp.plot

resp.plot2 <- lat_resp_dat2 |> 
  group_by(region, dev.temp) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~region); resp.plot2

resp.plot3 <- lat_resp_dat2 |> 
  group_by(region, dev.temp) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic(); resp.plot3

03. CTmax_analysis_resp_analysis

Elliott Schmidt

2024-05-23

Scenario

Load packages

Set working directory

Import data

Data manipulation

Exploratory data analysis

2nd order polynomial

3rd order polynomial

Fit model [random factors]

Priors

Models

Null

Model1

Model2

Model3

Model4

LOO

Fit model [fixed factors - polynomials]

Prior investigation

Priors

Models

Model1.p1

Model1.p1

DHARMa

MCMC diagnostics

model1.p2

model1.p2

DHARMa

MCMC diagnostics

model1.p3

model1.p3

DHARMa

MCMC diagnostics

LOO

Fit model [fixed factors - GAM]

GAM model k=4

GAM model k=4

DHARMa

MCMC diagnostics

GAM model k=3

GAM model k=3

DHARMa

MCMC diagnostics

Partial effects plots

conditional_effects

Model investigation

summary

R2

tidyMCMC

gather draws

bayesplot

emmeans - pariwise

emtrenns - pariwise

probabilities

summary figure